Hexagonal, square and stripe patterns of the ion channel density in biomembranes 
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Transmembrane ion flow througli channel proteins undergoing density fluctuations may cause 
lateral gradients of the electrical potential across the membrane giving rise to electrophoresis of 
charged channels. A model for the dynamics of the channel density and the voltage drop across the 
membrane (cable equation) coupled to a binding-release reaction with the ceU skeleton (P. Fromherz 
and W. Zimmerman, Phys. Rev. E 51, R1659 (1995)) is analyzed in one and two spatial dimensions. 
Due to the binding release reaction spatially periodic modulations of the channel density with a finite 
wave number are favored at the onset of pattern formation, whereby the wave number decreases with 
the kinetic rate of the binding-release reaction. In a two-dimensional extended membrane hexagonal 
modulations of the ion channel density are preferred in a large range of parameters. The stability 
diagrams of the periodic patterns near threshold are calculated and in addition the equations of 
motion in the limit of a slow binding-release kinetics are derived. 

PACS numbers: 89.75.Kd,87.16.Uv,05.65.-|-b,47.20.Ky 



I. INTRODUCTION 

Spatio-temporal pattern formation is ubiquitous in sys- 
tems driven away from thermal equilibrium [E S S 13 ■ 
Many physical, chemical and biological systems display 
dissipative structures, even though the underlying pat- 
tern forming mechanisms are often completely different. 
Nevertheless many of these patterns, especially those 
emerging at the primary bifurcation, belong to a few uni- 
versality classes [l| and patterns occurring in rather dis- 
parate systems share qualitative and unifying properties. 

Pattern forming processes in biological systems such 
as the fluid mosaic model, dilute filament-motor solu- 
tions (see e.g. 0, Iffl Si); actively polymerizing fila- 
ments [3, spiral waves in the cardiac system [13, skin 
patterning of the angle fis h nM or oscillatory dynam- 
ics in cell division [13, Il3l ll4L Il5| , are in general more 
elaborate than in classical pattern forming systems as 
for instance in fluid dynamics Q. In the latter case the 
equations of motion can be derived by using elementary 
conservation laws and phenomenological transport laws 
and accordingly, for various patterns in fluid dynamical 
systems a quantitative understanding has been achieved 
with a high precision [l|,|3|- These achievements can serve 
as a guide for the analysis of more involved biological or 
chemical pattern forming systems, where the respective 
models cover only the key steps of the complex biochem- 
ical reaction cycles. 

In the present work we investigate pattern formation of 
ion channels embedded in a biomembrane. Membranes 
are an important building block of living cells and play 
a key role for the biological architecture. They consist 
of a lipid bilayer which is rather impermeable and build 
the barrier to the cell environment. All the vital compo- 
nents needed inside the cell are transported across mem- 
branes through specific proteins. Especially for the sig- 
nal distribution along nerve cells (axons) the transport 
of ions through ion channels embedded in the membrane 
is essential since the transmembrane ion conductance is 



governed substantially by these discrete channels. The 
channel proteins are considered to move freely along the 
fluid lipid bilayer which is referred to as the fluid mosaic 
model 16], a concept that has attracted considerable at- 
tention. Describing the dynamics of ion channels within 
this framework, one finds transitions to various station- 
ary as well as time-dependent density patterns with pos- 
sible biological implications ;17. 18. 19. 20. 21. 22. 23. 24j. 
The binding-release reaction removes the conservation of 
mobile ion channels and as a consequence causes pattern 
forming instabilities with a finite wave number. Accord- 
ingly one expects in two-dimensional extended systems 
beyond a stationary bifurcation either stripes, squares 
or hexagonal patterns as prototype patterns. Here we 
focus on the competition between these patterns in the 
presented model system. 

Besides the fluid-mosaic model the channel concept 
[23 is the second central physical idea in the field of 
biomembranes. It was found that including their electro- 
diffusion properties |2a . Wx |2^ ion channels have an 
intrinsic propensity for self-organization [l7|- When a 
concentration gradient of salt across the membrane ex- 
ceeds a certain threshold, the conserved number of freely 
movable ion channels may organize into transient pe- 
rio dic p atterns which finally decay into global clusters 

The model of a fluid-mosaic of ion channels is only ele- 
mentary and neglects at least three important properties 
of real biomembranes: (A) An interaction with signal 
molecules may induce a reversible molecular transition 
which opens or closes an ion channel J25l | , (B) an interac- 
tion with the cell skeleton may immobilize ion channels 
|29| and (C) the excluded volume interaction between the 
ion channel molecules. The spatially dependent mobility 
of ion channels due to rafts [30| is also an effect neglected 
in the fluid mosaic model, but this heterogeneity effect is 
beyond the scope of the present work. 

The opening-closing reaction keeps the number of mo- 
bile ion channels conserved and its effects on the insta- 



bility of the homogeneous ion channel distribution have 
been investigated thoroughly in two recent publications 
[22, |2J| . Since membrane deformations coupled to the 
underlying cell skeleton {actin cortex) may also open or 
close ion channels ,31^ ^J , additionally we take here both 
the immobilization and the closing of the channels into 
account J22|. For the sake of simplicity we choose a 
model which combines the two processes: We consider 
a reversible binding-release reaction of ion channels with 
the cell skeleton and assume that this interaction induces 
a closing of the channels. Alternative models for pat- 
tern formation along the cell membrane take additional 
intermediate ste ps o f the opening-closing dynamics into 
account [l^ |22, 123, 1^ but a thorough analysis of the 
pattern formation processes in two spatial dimensions is 
not available yet. 

In the model we propose the closed ion channels which 
are considered to be bound to the cell skeleton are acting 
as a source for the mobile and open channels and there- 
fore, the free and open ion channels are not conserved 
anymore in contrast to previously discussed models. We 
find in this model different kinds of selforganization of the 
mobile ion channels: With a considerable binding-release 
reaction one has an unconserved number of open chan- 
nels and one finds (a) stable stripe or hexagonal patterns 
above threshold. This formation of stationary periodic 
patterns belongs to the same universality class as, for 
example, convection rolls in hydrodynamic systems, (b) 
The transition into the periodic pattern is either sub- or 
supercritical depending on the equilibrium constant, the 
relaxation time of the binding-release reaction and also 
on the strength of the excluded volume interaction of 
the ion channel molecules. In the limit of small binding- 
release reaction rates the model shows a crossover be- 
tween pattern formation for an unconserved and a con- 
served order parameter. For this crossover regime a re- 
duced equation is derived in Sec. IIVCI 

This work is organized as follows: In Sec^jwe describe 
the model system and we give the basic equations for the 
analysis in the subsequent sections. The linear stability 
of the homogeneous distribution of the ion channel den- 
sity and the onset of the patterns is discussed in Sec. lIIII 
The amplitude equations, describing the weakly nonlin- 
ear behavior of stripes, hexagonal and square patterns 
are derived in Sec. IIVI In this section also the nonlinear 
competition between these patterns is investigated by a 
thorough analysis. In Sec. numerical solutions of the 
model equations are presented. Those numerical results 
provide an estimate for the validity range of the pertur- 
bational analysis given in Sec. IIVI Concluding remarks 
and a discussion of the results are provided by Sec. I VII 



in close contact to another cell or to a membrane ca- 
ble as it occurs in dendrites and axons of neurons. In 
the first case the thin layer is given by the extracel- 
lular cleft and the bulk by the cytoplasm. A particu- 
larly important biological example of this case is the post 
synaptic membrane of a neuronal synapse. In the second 
case the narrow cylindrical cytoplasm plays the role of a 
one-dimensional cleft opposed by the extracellular bulk 
medium. 




FIG. 1: A fluid membrane that separates a narrow cleft 
of electrolyte from a bulk electrolytic phase is considered. 
Membrane proteins are embedded in the lipid bilayer. They 
are mobile along the membrane (diffusion coefficient D), they 
form selective ion channels across the membrane (conductance 
A), they bear an electrophoretic charge q and they interact 
with a filamentous substrate of the membrane (cell skeleton) 
via a binding-release reaction with rate constants Ko and Kc- 
Binding closes the channels by a conformational change. The 
system is driven by a concentration gradient of those ions 
which are conducted by the channels (Nernst-type potential 
E). The model refers to the biological situations of a cell- 
cell-contact (post-synaptic membrane of a synapse) and of a 
cylindric cellular cable (neuron dendrites). The cleft of the 
model corresponds to the extracellular space in the first case, 
whereas it corresponds to the narrow cytoplasm in the second 
case. The structure of the model is described by the density of 
mobile channels n(r, t) and by the voltage in the cleft v{r,t) 
as a function of space r and time t. 

Here we assume that the ion channels interact with the 
underlying cell skeleton, cf. Fig. [2 and Refs. [Mill Hi- 
In a reversible binding-release reaction among the ion 
channels and the cell skeleton the ion channels undergo 
also a conformational change and switch between an 
opened and a closed state. This binding-release reaction 
of the ion channels (IC) is described by a simple reac- 
tion scheme with two rate constants Kq and Kc, i.e. by 
an equilibrium constant Kbr — Kq/kc and a relaxation 
time tbb, = (kq + Kc)""^: 



Tf^bound ^ jr^free 

closed ^ open 



(1) 



II. MODEL SYSTEM 

We consider a model membrane with embedded ion 
channels separating a thin electrolytic layer from an elec- 
trolytic bulk medium. This may refer to a cell membrane 



The closed ion channels are bound to the cell skeleton 
and represent a source for free and open channels via 
the binding-release reaction. Accordingly the number of 
open and free channels is not conserved in our model, 
whereas in several previous investigations the number of 
free and open channels was conserved. Further mecha- 



nisms with similar consequences as the nonconservation 
of the ion channel number are also known [35J . 

The free ion channels with an electrical conductance A 
undergo a Brownian motion along the membrane with a 
diffusion coefficient D. They are selective for ions pre- 
senting a concentration gradient across the membrane 
which is described by a Nernst-type potential E. The 
proteins bear an effective electrophoretic charge q leading 
to a drift motion in a lateral electrical field. The current 
through the mobile and open channels and through a ho- 
mogeneous leak conductance of the membrane affects the 
local voltage in the cleft. In combination with an inho- 
mogeneous distribution of the ion channels this gives rise 
to lateral gradients of the voltage. 



A. Basic equations 

In a mean field approximation the local density n{r, t) 
of free and open channels (particles per unit area) is de- 
termined by diffusion and electrophoretic drift of the ion 
channels, the binding-release reaction and the local inter- 
action forces. The homogeneous density n is kept con- 
stant by a reservoir ric = const, of bound channels by 
an equilibrium binding-release reaction given by Eq. (Q 
with n = KsbjIc- The equations of motion for the chan- 
nel density may be expressed in terms of the deviation 

n — n ~n (2) 

from the mean density n, where the dynamics is com- 
posed of a lateral current gradient and a source-sink con- 
tribution 



n of the following from 



dth 



-V-j(n) 



tbr 



(3) 



The current depends on the part of the chemical poten- 
tial, which is independent of the electric field, and on the 
mobility of the charged channels times the local electric 
field 
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(7) 



The higher order contributions become important espe- 
cially for a negative diffusion constant Z), if a demixing 
between lipid proteins and ion channels in the membrane 
takes place. Since the diffusion constant is assumed to be 
always positive for the present problem, the higher order 
contributions may be neglected in most situations. How- 
ever, if the amplitudes of the spatial ion channel density 
modulations, as induced by the pattern forming mecha- 
nism discussed in this work, become strong, the second 
and third order terms in Eq. (jjj) , describing the effects of 
excluded volume interactions become important in some 
range of parameters. Here we investigate exemplarily 
the stabilizing effects of the fourth order contribution to 
Eq. (|ZI), i.e. with g2 ~ 0, ^ = but g ^ 0- In this case 
the equation of motion for n takes the following form 



dth =V^ (Dfi 



gn^) + 



qD_ 
kuT 



[{h + n) Vv] . 

Tbr 

(8) 



The voltage v in the cleft is obtained from Kirchhoff's 
law for each element of the membrane cable. Taking into 
account the current across the membrane and along the 
core of the cable or below a flat membrane, we obtain the 
Kelvin equation with the membrane capacitance C either 
per unit length in a one dimensional model or per unit 
area for a two dimensional model, respectively, the resis- 
tance R of the cleft per unit length (area) and the leak 
conductance G of the membrane per unit length (area) 

Cdtv^^\/'^v~Gv~Kn[v-E). (9) 

R 

The special case tbb. -^ oo and g = of these equations 
has been investigated in Refs. [i3, llM 123 • 



j(n) ~ — V/i — vqnVv 



with the mobility 



D 



(4) 



(5) 



the diffusion constant D, the effective charge q per chan- 
nel, and the local voltage drop across the membrane 
^(r, t). The charge q is assumed to be an effective charge 
taking into account screening as well as electro-osmotic 
effects J33 • The remaining part of the chemical potential 
may be derived from the free energy 



_5T_ 
Sh 



(6) 



where the free energy (per unit area S) of the channel- 
channel interaction is up to leading order in the deviation 



B. Scaled Equations 

We rescale Eq. ^ and Eq. ^ by introducing di- 
mensionless coordinates for space r' = r/A and time 
t' = t / T with the typical length scale of an electrical 
perturbation A = [R{An + G)]^^^^ and the time con- 
stant of displacement r = }? /D. We use normalized 
variables for the particle density N = {n — ri)/n and volt- 
age V — {v — vji)q/kBT with the resting voltage vb, = aE 
and the density parameter a = An/(An + G). Introduc- 
ing the normalized relaxation time Ty = RGD we obtain 
the normalized reactive Smoluchowski-Kelvin equations 

M 

dt.N =V'2 [N + gN^] + V' • [(1 + N)\I'V] ~ f3N , 

(10a) 

Tvdt'V = [V'2 -l]V -a{l- a) eN - aNV . (10b) 



The dynamics of the system is controlled by the fol- 
lowing three parameters: (i) The density parameter a 
characterizes the equilibrium of the binding-release re- 
action; (ii) the rate parameter f3 = t/tbr which char- 
acterizes the dynamics of the binding-release and the 
simultaneous opening-closing reaction; (iii) the control 
parameter e = —qE/{kBT) which characterizes the dis- 
tance to thermal equilibrium. Since the spread of the 
voltage is fast compared to the diffusion of ion channels 
(R^ 10*rj,C « liMF/cm^^D « Q.lnm^/s ^ r„ < 1 
J39j) we put Ty = in the following. For simplicity the 
primes of the new coordinates r' and t' are suppressed 
further on. 



III. THE ONSET OF PERIODIC PATTERNS 

The onset of spatial patterns takes place in a parame- 
ter range, where the homogeneous density n = n of the 
mobile channels, i.e. A^ = and V = Q, becomes linearly 
unstable with respect to small inhomogeneous perturba- 
tions. In order to calculate this instability the linear part 
of Eqs. (|10|) is transformed by an ansatz 



6N 



(Tt-\-ik.r 



(11) 



into linear algebraic equations. The solubility condi- 
tion of these equations determines for finite perturbations 
5N, 5V ^ and for Ty — Q the dispersion relation a{k) 



a(l 



a]e 



1 + P 



(12) 



with k — |k|, which is shown for three different values 
of the control parameter e in Fig. [^Ja). A perturbation 
grows in the range of the wave number k where the real 
quantity (j{k) becomes positive. The neutral stability 
condition a{k)—0 applied to the expression in Eq. 112|) 
gives the neutral curve as follows 




FIG. 2: Part (a) shows the dispersion relation u{k) as given 
by Eq. H12|l for a supercritical, a critical and subcritical value 
of tire reduced control parameter rj and for the rate parameter 
l3 = 0.08. Part (b) shows for four different values of the rate 
parameter /3 = 0.01,0.1,0.5, 1.0 (bottom to top) the neutral 
curves q(1 — a)eo(fe) (solid lines) with £o(fc) given by Eq. 1131 . 
The dashed curve in (b) marks the location of the minimum of 
these neutral curves, i.e. the critical wave number {kc,ec) (/3) 
as given by Eq. 11411 and Eq. 11151 . 



e (k) =. ^^ (l A 

a{l — a) \ fc2 



(13) 



where £o(fc) separates the range of stable from the unsta- 
ble parameter values. A set of neutral curves £o(fc; a, P) 
is shown in Fig. E^b) for different values of the rate pa- 
rameter (3. 

The minimum of eQ{k]a,f3) defines the critical wave 
number 



and the critical control parameter 
_ (1 + v^)' 



a(l — a) 
at which the basic state becomes first unstable. 



(14) 



(15) 



For £-values above the neutral curve eQ{k,a,(i) the 
growth rate a is positive and it takes its maximum at 
the wave number km 



ki = ^i + {i + ^) yrr 



whereby the reduced control parameter rj 



(16) 



(17) 



has been introduced. For a vanishing rate parameter j3 = 
the number of channels is conserved and the dispersion 
in Eq. I|12|l is at small values of k proportional to k^ which 
is similar to hydrodynamic excitation modes. To some 
extent this limit has already been investigated previously 

miiiiii. 



IV. WEAKLY NONLINEAR ANALYSIS 

For a finite value of the binding-release reaction param- 
eter (3 the fraction of the open ion channels is not con- 
served. In addition the homogeneous distribution of the 
ion channel density as well as the homogeneous voltage 
drop across the membrane may become simultaneously 
unstable against perturbations above a certain thresh- 
old. The perturbation with the wave number kc — /3^'^ 
has the largest growth rate, as described in the previous 
section. 

Near threshold and in two spatial dimensions this pe- 
riodic instability may lead to stripe, square and, if the 
up-down symmetry is broken, also to hexagonal pat- 
terns [y. In a parameter range where the amplitudes 
of these patterns are still small, the slow spatial vari- 
ations of these patterns may be described in terms of 
generic amplitude equations, a method used for many 
other physi cal, chemical and biological pattern forming 
systems [l|, 14^ ■ These generic equations and their cor- 
responding functionals are derived in this section from 
the present model system whereby the detailed scheme 
of derivation is given exemplarily for stripes in App. ^ 
The parameter ranges where each pattern realizes the 
lowest functional value and where two patterns coexist 
are determined in Sec. IIV IJI Far beyond the threshold 
the solutions are determined numerically and the ques- 
tion of preference of patterns is addressed by numerical 
simulations in Sec. 

In the range where the binding-release reaction be- 
comes rather slow and the number of channels nearly 
conserved, i.e. (3 (x rf, another set of equations is derived 
and discussed in part IIV CI The analytical solutions de- 
rived in the two cases /3 oc 0(1) and (3 o^ rf are compared 
with numerical solutions of Eqs. H1(J|) in Sec. 



A. Periodic patterns for a finite binding-release 
reaction, j3 oc 0(1) 

1. Stripe patterns 

For finite values of (3 the solution of the linear part of 
Eqs. 1)10(1 is in the simplest case spatially periodic in one 
direction. The two fields N and V may be written as a 
vector u = (iV, V) in the form 



uo =eoAe 



'kcT 



with the eigenvector 



eo = 



1 



En = -[l 



(18) 



(19) 



where c.c. denotes the complex conjugate of the preced- 
ing term. Due to rotational symmetry the orientation of 
the wave vector 



Kc — r^c 



(20) 



may be chosen parallel to the x-direction. Spatial vari- 
ations of the pattern which are slow on the length scale 
27r/fcc, are described in Eq. ((18|l by a spatially depen- 
dent amplitude A{x, y, t). The evolution equation for this 
spatially and time dependent amplitude is a so-called am- 
plitude equation, namely the Ginzburg-Landau equation 
in our case, which takes for an isotropic two-dimensional 
system, as for instance for the two-dimensional flat mem- 
brane, the following universal form pll3.l4(]||: 



TodtA 



v+eo[d.-,^d^y 






l\Af 



A. (21) 



The values of the coefficients To, £oand 7 depend on the 
specific pattern forming system 0,0, Eil. The relaxation 
time To and the coherence length ^0 in the amplitude 
equation H21|l may be derived from the dispersion rela- 
tion o'(fc) in Eq. H12I) and from the neutral curve eo{k) 
specified in Eq. H13|l by a Taylor expansion of both for- 
mulas around the critical point (fcc^c) (see for instance 
Refs. Qilil): 



da 
de 



^2 - 



1 d^eo 
2£c dk^ 



(22) 



Using the expressions in Eq. H12|) and Eq. \i'6\ we ob- 
tain their explicit form for the present model: 



To 



V^(l + V^) 



a 



{l + VP)' 



(23) 



The sign of the nonlinear coefficient 7 determines, 
whether the transition to the periodic state as described 
by Eq. (fTH|l is supercritical (7 > 0) or subcritical (7 < 0), 
and it may be derived by a perturbation calculation from 
the basic equations H10|) : 



7 = 



3.9 



1 



l + VP 



6a' 



{2 + 2^-ay 



1 + V^ 



(24) 



3^/? 



(4V/?-2a + l) (V^- 2a -1-1 



The common scheme for the derivation of 7 may be found 
for instance in Refs. [ll liol |4^ and the details of this 
calculation for the present system are given in App. ^ 

The linear coefficients in Eq. (|23|l depend only on the 
rate parameter /3. By increasing the binding-release reac- 
tions the relaxation time and also the correlation length 
becomes smaller. The nonlinear coefficient depends be- 
sides the rate parameter (3 also on the density parameter 
a and on the nonlinear interaction parameter g. In the 
limit P ~f for a ^ 1/2 the relaxation time To and the 
nonlinear coefficient 7 diverge. This behavior reflects the 
fact that in this limit the validity range of the amplitude 
equation is left and a different perturbation expansion 
has to be used as described in Sec. IIV CI below. 



The amplitude equation (|21|l for a stripe solution can 
also be derived from a functional T 



T^dtA 



5Ts 
5 A* 



(25) 



[A* denotes the complex conjugate of A) of the following 
form 



^'4/,* 



l\Af-lW + il 



S.-^fi\A 



(26) 



In the following we will focus on spatially homogeneous 
patterns and their competition, i.e. A{x^€) — > A{£). So 
the Ginzburg-Landau equation H21fl reduces to a Landau 
equation 



T^^tA = ^A-^\A\^A 
with a simplified functional 



(27) 



^^4//' 



7 



7i 



L\At-r^\A\^\=l^\At-r^\A\^. (28) 



Besides the trivial solution yl = 0, Eq. H27|l has a second 
stationary solution 



A^ 



(29) 



This solution exists in the supercritical case 7 > only 
above the threshold and for 7 < only in the range 77 < 
on the unstable branch of the subcritical bifurcation. As 
the herein presented expansion breaks down for subcriti- 
cally bifurcating stripes, higher order terms with respect 
to A would have to be taken into account in order to 
achieve a limitation of the amplitude A which may how- 
ever be determined by solving Eqs. (|10l) beyond threshold 
numerically as done in Sec. 

The parameter range of the supercritical and sub- 
critical bifurcation are separated by the tricritical line 
7(a,/3) = 0. This line is shown in Fig. |31in the a — fi 
plane for different values of the excluded volume param- 
eter (7, where the supercritical range may be extended by 
increasing the nonlinear parameter g. 

In the limit /3 — + and a -^ 1/2 the nonlinear coeffi- 
cient 7 = 3g+l/4 > is positive and the bifurcation still 
supercritical. Otherwise 7 diverges in the limit P —> 0. 
This is also in agreement with an alternative perturba- 
tion analysis as described in Sec. IIVCI 



2. Square patterns 

Square patterns can be described by a superposition of 
two periodic waves as given by Eq. ((TH|l 



uo = eo [Aie 



(kir 



^26^''^'") 



(30) 




FIG. 3: The lines describe the tricritical point, i.e. 
^{a,f3,g) = 0, for the stripe pattern in the a — /3 plane and 
for different values of the interaction coefficient g, respectively. 
On the right hand side of each curve corresponding to different 
values of g, 7 is positive and stripes bifurcate supercritically. 



whereby the two wave vectors ki and k2 have the same 
length but both are orthogonal to each other 



ki = kc 



and k9 = K- 



(31) 



Similar as for stripes, one may derive by a perturbation 
calculation from the basic Eqs. UlUI) the following two 
coupled equations for the two amplitudes Ai and A2 



TodtAi ^T]Ai- (7 I All 



x|A2p)Ai, (32a) 

TodtA2^vA2-{l\A2\^ + x\Ai\^)A2, (32b) 



wherein the spatial dependence of the amplitudes has al- 
ready been discarded. Herein tq and 7 are defined by the 
same expressions as for stripes in Eq. (|23|1 and Eq. H24|l . 
For squares one obtains with a perturbation calculation, 
similar as described in App. ^ the same expression for 
7 as in Eq. H24|) and the nonlinear coupling term x 



X = - 4 (2(1 - 2a) + V^j 

2(3g + Q)^~2(l-2a)^) 4(1 - 2a)^ 

^ 1 + vp 7m+W) 



(33) 
(34) 



As for stripes, the two coupled equations may be once 
again be derived from a functional 



TodtAi = - 



SA* 



(35) 



by determining the extremal value of the functional 

2 
•^0 = E (i l^^l' - ^ l^^l') + x\Ai\' IA2P . (36) 



i=l 



Apart from the trivial solution Ai = A2 = 0, the coupled 
amplitude equations H32|) have two types of stationary 



solutions of finite amplitudes. The first type corresponds 
to simple stripe solutions with only one finite modulus 

l^il- J^,|A2|=0 or \A,\^0,\A2\^J^ . 

(37) 

For the second type of solutions the amplitudes have 
identical moduli 



Ui 



\A2 



V 



7 + X 



(38) 



which corresponds to a square pattern as can be seen 
fromEq. ^ . 

If the sum of the two nonlinear coefficients is positive, 
i.e. 7 -|- X > 0, a square pattern bifurcates supercritically 
from the homogeneous state. A vanishing sum 7 + x = 
marks the tricritical line of the square pattern, the bifur- 
cation changes from a supercritical to a subcritical one. 
The tricritical line is displayed for different values of the 
interaction parameter g in Fig.^ where increasing values 
of g broaden also the range of supercritically bifurcating 
squares in the a — /? plane. 




FIG. 4: The nonlinear coefficient 7 + x as given by Eq. H33|l 
is positive and square pattern bifurcate supercritically in the 
range enclosed by the solid line for g — 0.0, for g = 0.1 in the 
range enclosed by the dashed line and for g — 0.5 between the 
dotted line. 



3. Hexagonal patterns 

In two-dimensional systems close to the threshold and 
without an up-down symmetry N,V — > —N,—V as in 
the Eqs. pOfl hexagonal structures are often preferred in 
some parameter range to stripe or square patterns [ij. 
In this section the amplitude equations of hexagons are 
presented as they result by a perturbation calculation 
from Eqs. H1U|) similar as described in App.^for stripes. 

Close to threshold a hexagonal pattern can be de- 
scribed by a superposition of three plane waves (stripe 



patterns) as given by Eq. H18(l , but where the three wave 
vectors enclose an angle of 27r/3 with respect to each 
other. The solution may thus be represented by 



uo = eo {Aie 



ikir 



A2e'^^'' 



Ase'"'"') + • 



(39) 



whereas the wave vectors ki (i = 1, 2, 3) are given by 



ki = kc 



and \i.2,3 = — 



-1 

±V3 



(40) 



The coupled amplitude equations for the three envelope 
functions Ai {i — 1,2, 3) are of the following form 



TodtAi =ijAi+SA*Al 



(^\A,\' + p\A2\^ + p\A,f)A, 



TodtA2 ='nA2 + SA*3Al 

-(7|A2P+H^3P 

TodtA3^T]A3 + 6AlA; 

-(7|A3p-fH^iP 



-p\AinA2 



-P\A: 



n^3 



(41a) 
(41b) 
(41c) 



with A* being the complex conjugate of Ai. tq and 7 are 
defined by the same expressions as for stripes and squares 
above in Eq. (|23|) and Eq. H24|l and the two nonlinear 
coupling constants d and p read within the scope of our 
model system 



(5 = 



P = 



These three coupled nonlinear equations H41|l nray again 
be derived via the relation 



l + y/j3-2a 

1 + VJ3 ' 
6.g 4a2 


2a{l + a) 

'^ {i + vpy 

2a{l - 2a) 


3(2 H 
4(H 

-|-3a. 


(42a) 


1 + v^ {1 + vpr 

(l-2a)(3 + 2a) 
4V^(1 + VP) 


(42b) 



TodtA, 



5T, 



H 



5A* 



(43) 



from a functional 
3 



^--EUi^ 



i=l 



'-V\A, 
5 {A1A2A3 



I A- 



AlA^Al) . (44) 



Eqs. (|41|l admit two types of homogeneous solutions. The 
first one corresponds to a stripe solution with only one 
non-vanishing amplitude. For hexagonal solutions the 
moduli of the three amplitudes |Ai| = |^2| = |^3| = A 
coincide, but if one allows still a relative phase shift (pi 
{i = 1,2,3), with 

A=Ae''^\ (45) 

one obtains the nonlinear equation 

= 7]A + 5 e"^A^ - (7 + 2p)A^ , (46) 



8 



with the sum of the three phase angles <I> = (/)i 
There are two real solutions of Eq. (|46|l 



A+ = 



1 



2(7 + 2p) 



5±^5^+ A-q (7 + 2p) 



(47) 



with A+ corresponding to the larger amplitude for 5 > Q 
and A- for ^ < 0. For 5 > the phase angle is $ = 0, 
which corresponds to regular hexagons and for 5 < the 
angle is $ = tt, which corresponds to inverse hexagons. 
Comparing for both solutions the functional T given by 
Eq. 144() . one finds that regular hexagons have the lower 
functional, i.e. T'^ < !F^, in the range with S > and 
inverse hexagons in the range oi S < with T^ > T^. 

The bifurcation from the homogeneous distribution of 
ion channels to a hexagonal modulation of the channel 
density is subcritical according to the quadratic nonlin- 
earity A^ in Eq. H46|) , which originates from the quadratic 
nonlinearity A* A* in Eqs. H41|) . However, the amplitudes 
Ai are still bounded by cubic nonlinearities in the param- 
eter range of a positive nonlinear coefficient 7 + 2/3 > 
in Eq. 146|l . This nonlinear coefficient vanishes along the 
lines shown for different values of the parameter g in the 
a — /3 plane in Fig.[Sl If this coefficient becomes negative, 
i.e. 7 + 2p < 0, Eqs. H41() do not have any stationary, 
finite amplitude solutions. In this case one needs either 
a higher order expansion or the amplitudes of hexagons 
have to be determined by solving the basic equations H1U|) 
numerically. Increasing values of the nonlinear interac- 
tion parameter g enlarges the parameter range wherein 
stationary solutions of the form (|47|l occur. 




FIG. 5: The nonlinear coefficient 7 -|- 2p for hexagons is pos- 
itive between the solid line for g = 0.0, for p = 0.1 between 
the dashed line and for g = 0.5 between the dotted line. In 
each range the amplitude of the hexagonal solution is limited 
by the cubic terms in Eq. I|41^ . 



B. 



Competition between patterns 



hexagons is simultaneously limited by a cubic term. This 
range becomes even larger with increasing values of g as 
indicated by the ranges in Fig. 13 Fig.^and Fig.|Sl So the 
interesting question arises, which of the three solutions 
is preferred in this overlapping range. 




In the shaded subrange in Fig. El stripes and squares 
bifurcate both supercritically and the amplitude of the 



FIG. 6: The bifurcation behavior of stripes, squares and 
hexagons is shown in the a-j3 plane for g — 0.5: In the shaded 
range the amplitudes for stripes, squares and hexagons are 
limited by a cubic nonlinearity, i.e. 7>0, 7-|-x>0 and 
7 + 2p > 0. On the right hand side of the dashed line, which 
is determined by 7 = 0, stripes bifurcate supercritically and 
squares do so in the range inclosed by the solid line, which is 
defined by 7 -I- x = 0. Between the dashed-dotted line which 
is determined by 7 -|- 2p = hexagons are limited by cubic 
order terms. 

One criterion is the comparison of the values of the 
functionals, i.e. to which solution belongs the lowest 
value of the respective functional J-'. The second criterion 
is the linear stability of each of the nonlinear solutions, 
i.e. in which subrange of the overlap range of parameters 
becomes one of the solutions linear unstable with respect 
to small perturbations. 



1. Comparison of the functionals for stripes, squares and 
hexagons. 

For one set of parameters, /3 — 0.15, rj = 0.06 and 
g = 0.5, the functionals of the three patterns are shown 
in Fig. as a function of a in the range where each of 
them bifurcates supercritically. This figure indicates, in 
which region the respective pattern has the lowest value 
of the functional. 

For small {a < 0.49) and very large values of a 
{a > 0.98) the functional !Fq has the lowest value, 
i.e. squares have the lowest energy and are accord- 
ingly preferred. Regular hexagons J^^ are preferred in 
the range 0.49 < a < 0.56, stripes J^s in the range 
0.56 < a < 0.86 and finally inverse hexagons T]^ in the 
range 0.86 < a < 0.98. These respective ranges change 
as a function of 77, (3 and g. 




expressions 



-0.003 



-0.004 



FIG. 7: The functional T per unit area is shown for stripes 
(jFs), squares (Tq), regular (T^) and inverse hexagons (T^) 
as a function of the parameter a as well as for the parameter 
values /3 = 0.15, ri — 0.06 and g = 0.5. 



Plotting the crossing points of the curves in Fig. [7| as 
a function of the kinetic parameter /3 leads to a phase 
diagram as presented in Fig. |S1 for 77 — 0.06 and g — 
0.5. In this figure hexagons have a lower functional value 
than stripes beyond the upper dashed line and below 
the lower dashed line and a lower functional value than 
squares between the dotted lines. Taking the competition 
between squares and stripes into account too, stripes are 
preferred in the dark shaded range, squares in the bright 
shaded range and hexagons in the medium shaded range. 
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FIG. 8: For r\ — 0.6 and g — 0.5 the parameter ranges are 
shown where stripes (dark), hexagon (medium) and squares 
(bright) have the lowest value for the functional T. Along the 
dotted line one has J-r = ^q , along the dashed line Th ~ ^s , 
and along the dash-dotted line Ts = Tq. 

Inserting the solutions of stripes in Eq. (|29|l and 
squares in Eq. (|38|l into their functionals in Eq. H28|l and 
Eq. (|36|l the functionals can be reduced to very simple 



:fs 



27 



v 



T, 



Q 



1 

7 + X 



?7 



(48) 



Hence the comparison of the functionals for squares and 
stripes is independent of the reduced control parameter 



Tq>Ts ^ x>7>0. 



(49) 



Accordingly the curves, cf. the dash-dotted line in 
Fig. |S1 as calculated from the condition 7 = x of equal 
functional values, separate the regions where the func- 
tionals of stripes or squares have the lower values. 

However, a comparison with the functionals for regu- 
lar and inverse hexagons is not independent of rj. Since 
hexagons bifurcate subcritically their amplitude is al- 
ready finite and they have lower functional values at 
threshold 7/ = 0, i.e. ^h{v = 0) < J-'s = Tq — 0. 
Hexagons are therefore always preferred close to the 
threshold. Stripes and squares are always favored with 
respect to hexagons beyond some critical values rj > 
ris{a.,(3) and 77 > r]Q{a,f3) which are determined by the 
conditions J^h = ^s and J^h = ^q, respectively. Ac- 
cordingly, with decreasing values of 77 the range in the 
a — (3 plane increases where hexagons have the lowest 
functional. As indicated by Fig. |51 the ranges of stripes 
and squares become smaller and smaller. For small values 
of 77 square patterns are suppressed nearly completely. 

For the interesting limiting case, a — 1/2 and /3 ^ 0, 
one has i5 = and p — Gg + 1/2. Thus one expects in 
two spatial dimensions for /3 <C 1 and a = 1/2 as well as 
small values of 7; a clustering of ion channels to stripes. 
For finite values of f3 and 7; > stripes are preferred 
to hexagons in the neighborhood of the curve ^ = 0, cf. 
Fig. 1^ This region broadens with increasing values of 77, 
as indicated by Fig. |51 



2. Linear stability analysis 

A comparison of the values of the functionals for the 
respective patterns is one criterion for their basin of at- 
traction. A linear stability analysis of the patterns shows, 
as described subsequently, that two patterns can both co- 
exist in a region around the parameter range where their 
functionals agree. 

a. Stripes vs. hexagons. Stripes are still linearly sta- 
ble even if their free energy is already higher than that 
of the hexagons, as can be seen by investigating the lin- 
ear stability of stripes with respect to small amplitude 
perturbations 6Ai 



Ai 



6A, 



Ao = SAo 



A, = &4, 



(50) 



Linearizing Eqs. H41|) with respect to the small functions 
6Ai(t) and solving the resulting linear differential equa- 
tions, one obtains the stability boundary as described by 
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FIG. 9: Dependence of the phase diagram of square, hexagon 
and stripe solutions (bright, medium and dark gray) on the 
control parameter rj: (a)-(d) rj — 0.6,0.4,0.2,0.1, g = 0.5. 
For small values of rj the region for hexagon solutions in the 
a-/9-plane increases and other solutions, squares quite more 
than stripe patterns, are suppressed. 



the condition (see e.g. Ref. |4J|) 



(51) 



By a similar stability analysis of the hexagonal solutions 
one obtains the following stability boundary |4J| 



^7(7-/3) -(7-2/9)^2 <0. 



(52) 



While stripes have the lower free energy between the solid 
lines in Fig.|Sl stripes become only unstable at the dashed 
line and the hexagons at the dash-dotted line in Fig. |S1 

h. Stripes vs. squares. By a linear stability analy- 
sis of the stationary solutions given by Eq. l|77|l and by 
Eq. (|38ll one finds that stable stripes are preferred in the 
range x > 7 > of the nonlinear coefficients and squares 
in the parameter range |xl < 7 !4vi- These ranges coin- 
cide with the ranges where both patterns have their lower 
functional values (cf. Eqs. (|49|) ). Therefore stripe and 
square patterns do not coexist. 

c. Squares vs. hexagons. Numerical results in 
Sec. IV Bl show that the amplitude equation for squares 
is a good approximation only for very small values of 77. 
In this range square patterns are nearly completely sup- 
pressed by hexagonal pattern. Thus the quite compli- 
cated stability analysis of square patterns vs. hexagons 
(see e.g. ;46]) has been left out. 




FIG. 10: Range of coexistence of hexagons and stripes for 
77 = 0.1 and g = 0.5. In the grey regions hexagons and 
stripes arc both linearly stable with respect to small ampli- 
tude perturbations. Therefore both patterns coexist in a fi- 
nite range around the dashed lines where their functionals 
agree, J-s ~ Th. Along the dotted line the nonlinear coeffi- 
cient 5(a, /3) = vanishes. The solid line indicates the region 
where both structures are limited by a cubic nonlinearity in 
their amplitude equations (7 > and 7 -I- 2p > 0). 



C. "Nearly conserved case" - (3 oarj 

The amplitude equations for stripes, squares and 
hexagons were derived in Sec. IIV bl under the assumption 
of a finite wave number kc — f3^^^ oc 0(1) and k^ 3> rj. 
In the limit of a conserved number of open ion channels, 
i.e. (3 = 0, beyond the onset of patter formation a clus- 
tering of channels takes place [43 , that may be described 
with a single component Cahn-Hilliard equation. In this 
section we show that in the limit of small values of the 
rate parameter /3 with fi o^ rf a. reduction of Eqs. (|10|l to 
a single equation is still possible, whereby the resulting 
equation covers the qualitative properties of demixing as 
well as of periodic pattern formation. 

From the previous section it is known that in the limit 
/3 -^ the bifurcation is subcritical in the range of small 
values of /3, besides a = 1/2. Therefore the appropriate 
expansion, which leads also to finite amplitude solutions, 
is the combined expansion of /3 ex 77^ and of a = 1/2 + 
a^/rj (with a ex 0(1)), as described in this section. 

a. Expansion in the regime (3 Oi rf' and a arbitrary. 
In this case we expand the two fields N and V with re- 
spect to powers of 77 in the following manner 



N =r]No + ri^Ni + 
V =r^Vn + ri^Vi + . 



(53a) 
(53b) 



In addition we also introduce slow space variables, 
X,Y = ^x,y/rjy, and a slow time variable T = rf't. 
Choosing [3 = rj^P, with /3 oc 0(1) and expanding 
Eqs. H1U|) with respect to powers of 77 one obtains a single 
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order parameter equation for A^o^ 



drNo 



V2 



mo 



V^ + l + 2V/3MVo 



" - 2 ) ^0 



(54) 



wherein V = {dxtdy) has been used. The voltage follows 
immediately via the identity 



Vn=-Nn 



(55) 



By rescaling space, time as well as the amplitude {N 
rjNo) this equation takes the form 



dtN 



V 



V2 + r/ + 2^^) V^TV, + (a - ^ ) ^(iV^) 



m, 

N . 



(56a) 
(56b) 



This equation for TV shares similarities with the damped 
Kuramoto-Sivashinsky equation 0,|43- The only differ- 
ence is in the nonlinearity y^{N'^) = 2{{VNy+NV^N) 
because the Kuramoto-Sivashinsky equation includes as 
nonlinearity (VTV)^ only. The additional term, NW'^N, 
however, changes the dynamics and stability of the solu- 
tions completely beyond threshold, 77 > 0. While one has 
for the Kuramoto-Sivashinsky equation "turbulent" but 
bounded solutions, the solutions of Eq. (|56a|l are always 
divergent according to the nonlinear diffusion NW'^N. 
For a — 1/2 the nonlinear coefficient in Eq. (|56a|) van- 
ishes completely, which suggests a different expansion 
close to this point, as described in the next paragraph. 

b. Expansion in the regime f3 (x rj"^ and (a — 2) oc 
y/rj. Here also the parameter a is expanded with respect 
to the small parameter r/: a = 1/2 + y/ffa with a ex 0(1). 
Expanding the fields iV and V with respect to powers of 



N =ri^No + T]Ni +ri^N2 + 
V ^ri^VQ + riVi +r]iV2 + .. 



(57a) 
(57b) 



yields from Eq. (|10|l at leading order a single equation 
for Nq, which after rescaling N — WffNo takes the form 



aiV = - V 



2 


(y^ + Tj + 2^^N- 


HV^ 


^+.)*' 


-m, 


(58a) 








(58b) 



V = -N 



The cubic nonlinearity now limits the amplitudes of the 
solutions to finite values, because 1/12+^ is positive even 
in the limit /3 ^ 0. For (3 — and a = 1/2, corresponding 
to a = 0, the transition to the inhomogeneous channel 
distribution is continuous while it is discontinuous for 
a 7^ but remains bounded. In the limit of a conserved 



channel density TV, i.e. 0^= 0, equation (|58a() is of the 
Cahn-HiUiard type [jgllsot . 

Eq. H58a|l covers all qualitative features of the basic 
equations IjlUI) in both cases, in the limit /? = and 
for /? ^ 0. Therefore, similar as in the previous section 
one can also derive the amplitude equations for stripes, 
squares and hexagons by starting from the modified equa- 
tion H58a|) instead of Eqs. (|10|) . These derivations are 
much simpler for Eq. H58a() but the results are now re- 
stricted to a range along the line a '^ 1/2, close to 
the threshold and to small values of the rate parame- 
ter /3 oc 77^. Thereby one obtains again the amplitude 
equation for stripes as given by Eq. H21() . for squares as 
by Eq. (|32|l and for hexagons as by Eq. (|41|) . but now 
with slightly modified expressions for the coefficients. 



To 



1 
1 



(5 = 1 



35, 



2a, 
1 



P = X = 



65. 



(59) 



These expressions may also be recovered from the formu- 
las given in Sec. lIVTO in the limit /^ — > and a -^ 1/2. 



V. NUMERICAL RESULTS 

The amplitude equations, as given for the present sys- 
tem in the previous sections, are exemplarily derived in 
the appendix by a perturbational calculation. However, 
the validity range of these equations and their solutions 
is a priori unknown. An estimation of this range can be 
provided by comparing the analytical solutions with nu- 
merical simulations of the basic equations (|10|l , as done 
in this section. 



A. Stripe patterns 

In the range of supercritically bifurcating stripe pat- 
terns determined analytically in the previous section, the 
analytical and numerical solutions of Eqs. (|1(J|I are com- 
pared in Fig. 111! as a function of the control parameter 
?7 for a = 0.4, /3 = 0.1, g = 0.5 and ry = 0.1. There is 
a fairly good agreement between the analytical and the 
numerical solutions up to about rj = 0.1. Since the stripe 
pattern is stationary it does not depend on the actual 
value of Ty used in simulations. 



B. Square patterns 

Numerical simulations do not show stable square pat- 
terns besides transients to finally hexagonal patterns. 
Choosing a special geometry with a very small system 
length (L = 27r/fcc) hexagonal patterns can be suppressed 
in numerical simulations. In this case the system shows 
the square patterns as predicted. By comparing the am- 
plitudes Ai as given for squares by Eq. H38|) for a = 0.4, 



12 



A 



0.20 



0.15 



0.10 



0.05 



0.00 




0.00 



0.02 



0.04 



0.06 



0.08 



1 



FIG. 11: The amplitude of a stripe pattern as determined by 
the amplitude equation (solid line) and by numerical solutions 
of Eq. HIO^ (circles) are compared as a function of the control 
parameter. For the parameters a = 0.4, 13 — 0.1, g = 0.5 
and TV = 0.1 the bifurcation was supercritical, i.e. 7 > 0. 



C. Hexagonal patterns 

Close to threshold hexagons are the preferred pattern 
in a wide range of parameters. In a range where stripes 
and squares bifurcate supercritically, but where hexagons 
are already preferred, at a = 0.4, /3 — Q.l and g — 0.5, 
the analytically, of. Eq. H47(l . and the numerically ob- 
tained solutions are compared in Fig. 1131 
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(3 = 0.1 and g = 0.5 with the numerical obtained solu- 
tion, as depicted in Fig. 1121 we find as a function of the 
reduced control parameter only an acceptable agreement 
for very small values of 77 < 10~^. But then we find 
that squares are only preferred according to the analyti- 
cal calculation in the range rj > 0.1, which is far beyond 
the validity range of the amplitude equations for squares. 
This is an explanation why we do not find the predicted 
squares by numerical solution of Eqs. H10|) . 
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FIG. 12: The amplitudes of squares, A\ = A2, are shown 
as a function of the reduced control parameter 77 and for the 
parameters a — 0.4, /3 = 0.1, g — 0.5. The solid line is 
given by Eq. 13811 and the circles are obtained from numerical 
simulations of the model equations. 



FIG. 13: The amplitudes Aj of the hexagonal pattern are 
given as a function of the reduced control parameter 77. The 
solid line corresponds to the analytical solution A+ and the 
dashed line belongs to the unstable solution A-, where both 
are given by Eq. 1471 . The data points are obtained from 
the numerical solution of Eqs. IIOII . The parameters (3 = 0.1, 
a = 0.4 and g — 0.5 have been used. 



D. Anharmonic solutions and clustering of ion 
channels for near conservation 



Increasing the control parameter up to i] = 0.9, far 
beyond the validity range of the amplitude equations, 
the density N{x) becomes rather anharmonic as shown in 
Fig-CUb). Closer to the threshold at 77 = 0.04 but in the 
range where stripes bifurcate subcritically and where the 
amplitudes take immediately large values, at a = 0.95 
and a = 0.03, the solutions are also very anharmonic as 
shown for N{x) in Fig. I14r a) and ll4f c). In both cases 
each peak in Fig. I14r a') and (c) takes already a similar 
shape that are typical for the clusters in the conserved 
limit /3 = 0. 



VI. DISCUSSION AND CONCLUSION 

In Ref. |22| a model for the dynamics of ion channels 
including electrophoresis, an opening-closing reaction as 
well as a simultaneous binding-release reaction has been 
introduced. Compared to this earlier work, we have taken 
into account the effect of the excluded volume interaction 
between ion channel molecules. In addition the analysis 
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FIG. 14: The stationary and normalized distribution of ion 
channels A''(a::) is shown for a stripe pattern above threshold 
for four different sets of parameters. In all cases /3 — 0.03 and 
<; = 0.1 were fixed and therefore also the critical wave number 
kc — 0.416. The system length is L = STr/kc- 



of the bifurcations beyond the threshold of pattern for- 
mation has been extended to two spatial dimensions. 

In terms of amplitude equations we give a detailed 
analysis of the competition between stripes, squares and 
hexagonal patterns. We have found that immediately 
above threshold hexagonal patterns are preferred in a 
large range of parameters, whereas further beyond the 
threshold stripes are preferred in an increasingly larger 
parameter range. 

The validity range of the amplitude expansion has been 
tested by solving the model equations numerically with a 
pseudo-spectral code. By this comparison we have shown 
that the amplitude equations for squares have a very 
small validity range. The expansion for squares breaks 
down two orders of magnitude earlier than the one for 
stripes and hexagonal structures. Accordingly, the large 
range where square patterns should be favored as pre- 
dicted by the amplitude expansion is not confirmed by 
the numerical analysis of the model equations while the 
amplitude equations describing the competition between 
stripes and hexagons are in a fairly good agreement with 
the numerical simulations close to the threshold. 

In the limiting case /3 — > implying a conserved num- 
ber of open ion channels the patterns have a strong sim- 
ilarity with those patterns occurring during ion channel 
clustering ji^l in systems with a binding-release reac- 
tion. Near threshold and in the limit of /9 ex 77^ the 
model equations can be reduced to a single model equa- 
tion which shares similarities with different variants of 
the Swift-Hohenberg equation [J as well as the Cahn- 
Hilliard equation [i^, IS^I . On the basis of this reduced 
equation essential effects related to the binding release 
reaction are already captured. 

It is an interesting question whether the curvature of 
membranes influences the pattern competition especially 
between stripes and hexagons. One expects that in such 



a system the effects of a broken up-down symmetry be- 
come stronger leading to an additional enlargement of 
the parameter range where hexagonal patterns occur. 

If the binding-release reaction is replaced by an 
opening-closing reaction, the formation of spatially pe- 
riodic either stationary or oscillatory patterns has been 
reported |23 . 12^ . There the nonlinear behavior of stripe 
patterns has been discussed only partially and only in 
one spatial dimension. A competition of patterns as de- 
scribed in the present work will be expected but also 
a competition between two-dimensional stationary and 
time-dependent patterns. 

Our calculations are related to experiments of the type 
as investigated in Ref. [SJ where ion channels are studied 
in vitro. The electro-diffusive model at hand seems to be 
relevant for in vivo systems as shown by experiments on 
the effect of electric fields on clustering of acetylcholine 
receptors [23, |23, 123 . In membranes composed of sev- 
eral types of lipids a selforganized structuring has also 
been described in terms of hpid rafts [Sfl. Therefore 
one expects in such cases a spatially varying mobility of 
proteins embedded in the membrane. How this hetero- 
geneity affects the formation of patterns is another inter- 
esting question as has been investigated for other model 
systems [Hill, HI HI, H^]. A detailed analysis of all 
these questions may be given in forthcoming works. 
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APPENDIX A: DERIVATION OF THE 

AMPLITUDE EQUATION FOR THE STRIPE 

PATTERN. 

1. Basic equations in Matrix- Notation 

Rewriting the basic Eqs. (|10|l to a matrix notation 

{M-dt + C)u = 'N{u), (Al) 

allows a more compact formulation of the derivation of 
the generic amplitude equations of stripe patterns as 
given by Eq. (|^ . of square patterns as given by Eqs. l|5^ 
or of hexagonal patterns as given by Eqs. (|41|l . The two 
components of the vector u are the normalized channel 
density N and the reduced voltage V, whereas the ma- 
trices Ai and C represent the linear parts of Eqs. H10|) 



M 



1 




C 



/3-V2 
a(l — a) 



-V2 
1-V2 



and the vector N the nonlinearity: 

N = 



V{NVV)+gV^{N^) 
-aNV 



(A2) 



(A3) 
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The linear coefBcients of Eq. H21() follow directly from 
the linear stability analysis as described in Sec. 11111 and 
Sec. IIV A II but the nonlinear coefficient 7 in the same 
equation is determined by the perturbation expansion 
described in the next subsection. 



Using the ansatz of Eq. (|A6|) the leading contribution 
Ni(uo) of N(u) has the explicit form 



Ni 



-So 






+ c.c. 



(A15) 



2. Nonlinear coefficient 

The sign of the nonlinear coefficient 7 in Eq. (|21() de- 
termines, whether one has a sub- or a supercritical bi- 
furcation to the periodic state given by Eq. H18|l . The 
scheme of the derivation of 7 and the respective ampli- 
tude equations may be found in various references, as for 
instance in Refs. 0, [i^, |43| • This scheme for the deriva- 
tion of 7 is summarized for the present system in this 
appendix. The starting point is an expansion of the so- 
lutions of Eqs. l|Aip with respect to powers of the reduced 
control parameter rj as given by Eq. H17|l: 



u(r,t) =r/2uo(r,t) +?7Ui(r, t) +r/2U2(r,t) -f 

Since the vector N is a nonlinear function of u. 
be also expanded with respect to powers of 77: 

N(u) =?7Ni(uo) + 77^N2(uo,ui) + ... . 



(A4) 
it may 

(A5) 



To the ansatz in Eq. H18|) for a homogeneous stripe pat- 
tern 



Uq = 60^6*'''='' +C.C 

a multiscale analysis in time is added 
dt—>dt + -qdT 



(A6) 



(A7) 



to account for the variation of the amplitude A — A{T) 
on a very slow time scale T = rjt. Together with the 
relation e = edl +ri) finally the basic equations in (|A1|) 
may be rearranged into a powers series with respect to 
772 leading to the following hierarchy of equations: 



(Mdt + £0) uo = , 
{Mdt + Co) u, = Ni{uo), 

{Mdt + Co) U2 = - (£2 + MOt) uo 



with the linear operators 
£0- 



A = 



a(l - a)ec, 



a{l — a)ec 



-V2 
1-V2 



and the nonlinearities 



Ni 



N2 = 



V (No^Vo) 
-aNoVo 

ViNoVVi + NiVVo) + gV^Nl^ 
-a{NoVi + NiVo) 



(A8) 
(A9) 

-N2(U0,Ui) 

(AlO) 



(All) 
(A12) 

(A13) 
(A14) 



The solution Ui of Eq. ljA9p has to be of the same form 
as the inhomogeneity Ni. This leads to the ansatz 

ui = (5ieo + B2ei)A^e^'^^'' + {BsSo + B^e^) \A\^ + c.c. 

(A16) 

using the two eigenvectors §0,1 = (1, £^0,1) of Co with 

i + VP 



Eo = -(i + Vp) , 



Ei^ 



VP 



(A17) 



After inserting of ljA16p in ljA9p one obtains by compar- 
ison of coefficients: 



n ^ vp{i + vp) v^ 



B2-U2^P-^,, 



B, 



i + VP'- 



Bi 



-B. 



(A18) 



At the next higher order, Eq. IJA10|) . one has to deal 
with the second order correction of the linear operator, 
£2 , and of the vector N2 as defined in Eq. (|A12(I and Eq. 
(IA14p . It is not necessary to solve Eq. IJA10|) explicitly. 
One can use Fredholm's alternative or one simply can 
take advantage of the following property 



(fne^- 



,£oU2) ^^ dr e-'''^''^CoU2{r,t) = 



of the left eigenvector 



fo 



^0 



Fo 



1 + V^' 



(A19) 



(A20) 



which spans the adjoint kernel of Cq. Since Uq and Ui 
have an explicit dependency only on the time scale T 
but not on t the corresponding derivatives can be ne- 
glected. Accordingly all the terms at the right hand side 
of Eq. IjAlOp projected onto fge*'^'^'' also have to vanish: 



(foe*=^ N2 - (£2 + XSt) Uo) = . 



(A21) 



This provides the solubility condition for the determina- 
tion of the amplitude A. For this purpose the contribu- 
tions to the expressions £2Uo and N2 which are propor- 
tional to e'^^"'' are collected. According to the Fredholm's 
alternative we obtain after projection 



-ylAl^A-A + rdrA^O 



(A22) 
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with the nonlinear coefRcient 
3.9 



7 



1 r 6^2 - (a - 2(1 + y^)y 

l + V/3^3 



1 + V^ 



+ — =(4v/^-2a+l)(V^-2a+l) 



and the relaxation time 



(A23) 



The nonlinear coefRcient depends on the rate parameter 
(3, the density parameter a and on the parameter g for 
the excluded volume interaction. The relaxation time of 
the pattern as well as the nonlinear coefficient 7 both 
diverge for a frozen binding- release reaction (/3 — > 0). 
In this limit the wave number kc tends to zero and the 
assumptions made for the derivation of the amplitude 
equation are not longer fulfilled. 



VP{i + VP)' 



(A24) 
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